%  
%  This program replicates the empirical work reported in 
%  "Testing Monotonicity of Mean Potential Outcomes in a Continuous Treatment 
%   with High-Dimensional Data" 
%

clear; tic;
rseed=RandStream('mt19937ar','Seed',10220408);
RandStream.setGlobalStream(rseed);

% Parameters
B=1000;
alphas=(0.001:0.001:1)';
eta=10^-6;
epsi=10^-6;
K_iter=5;
trim=0.0001;   % 0.0001
TXs=1;               % polynomials
K=5;
 
% Load data
for jj=1:64
    if jj==1; data0=readmatrix('jobcorpY1.csv'); c0=1; T_l=2000; T_u=3000; N=50; end  
    if jj==2; data0=readmatrix('jobcorpY2.csv'); c0=1; T_l=2000; T_u=3000; N=50; end      
    if jj==3; data0=readmatrix('jobcorpY3.csv'); c0=1; T_l=2000; T_u=3000; N=50; end  
    if jj==4; data0=readmatrix('jobcorpY12.csv'); c0=1; T_l=2000; T_u=3000; N=50; end              
    if jj==5; data0=readmatrix('jobcorpY1.csv'); c0=1; T_l=1000; T_u=2000; N=50; end  
    if jj==6; data0=readmatrix('jobcorpY2.csv'); c0=1; T_l=1000; T_u=2000; N=50; end      
    if jj==7; data0=readmatrix('jobcorpY3.csv'); c0=1; T_l=1000; T_u=2000; N=50; end  
    if jj==8; data0=readmatrix('jobcorpY12.csv'); c0=1; T_l=1000; T_u=2000; N=50; end      
    if jj==9; data0=readmatrix('jobcorpY1.csv'); c0=1; T_l=40; T_u=1000; N=50; end  
    if jj==10; data0=readmatrix('jobcorpY2.csv'); c0=1; T_l=40; T_u=1000; N=50; end      
    if jj==11; data0=readmatrix('jobcorpY3.csv'); c0=1; T_l=40; T_u=1000; N=50; end  
    if jj==12; data0=readmatrix('jobcorpY12.csv'); c0=1; T_l=40; T_u=1000; N=50; end  
    if jj==13; data0=readmatrix('jobcorpY1.csv'); c0=1; T_l=40; T_u=3000; N=50; end  
    if jj==14; data0=readmatrix('jobcorpY2.csv'); c0=1; T_l=40; T_u=3000; N=50; end      
    if jj==15; data0=readmatrix('jobcorpY3.csv'); c0=1; T_l=40; T_u=3000; N=50; end  
    if jj==16; data0=readmatrix('jobcorpY12.csv'); c0=1; T_l=40; T_u=3000; N=50; end      
    if jj==17; data0=readmatrix('jobcorpY1.csv'); c0=1; T_l=2000; T_u=3000; N=40; end  
    if jj==18; data0=readmatrix('jobcorpY2.csv'); c0=1; T_l=2000; T_u=3000; N=40; end      
    if jj==19; data0=readmatrix('jobcorpY3.csv'); c0=1; T_l=2000; T_u=3000; N=40; end  
    if jj==20; data0=readmatrix('jobcorpY12.csv'); c0=1; T_l=2000; T_u=3000; N=40; end              
    if jj==21; data0=readmatrix('jobcorpY1.csv'); c0=1; T_l=1000; T_u=2000; N=40; end  
    if jj==22; data0=readmatrix('jobcorpY2.csv'); c0=1; T_l=1000; T_u=2000; N=40; end      
    if jj==23; data0=readmatrix('jobcorpY3.csv'); c0=1; T_l=1000; T_u=2000; N=40; end  
    if jj==24; data0=readmatrix('jobcorpY12.csv'); c0=1; T_l=1000; T_u=2000; N=40; end      
    if jj==25; data0=readmatrix('jobcorpY1.csv'); c0=1; T_l=40; T_u=1000; N=40; end  
    if jj==26; data0=readmatrix('jobcorpY2.csv'); c0=1; T_l=40; T_u=1000; N=40; end      
    if jj==27; data0=readmatrix('jobcorpY3.csv'); c0=1; T_l=40; T_u=1000; N=40; end  
    if jj==28; data0=readmatrix('jobcorpY12.csv'); c0=1; T_l=40; T_u=1000; N=40; end  
    if jj==29; data0=readmatrix('jobcorpY1.csv'); c0=1; T_l=40; T_u=3000; N=40; end  
    if jj==30; data0=readmatrix('jobcorpY2.csv'); c0=1; T_l=40; T_u=3000; N=40; end      
    if jj==31; data0=readmatrix('jobcorpY3.csv'); c0=1; T_l=40; T_u=3000; N=40; end  
    if jj==32; data0=readmatrix('jobcorpY12.csv'); c0=1; T_l=40; T_u=3000; N=40; end         
   
    if jj==33; data0=readmatrix('jobcorpY1.csv'); c0=-1; T_l=2000; T_u=3000; N=50; end  
    if jj==34; data0=readmatrix('jobcorpY2.csv'); c0=-1; T_l=2000; T_u=3000; N=50; end      
    if jj==35; data0=readmatrix('jobcorpY3.csv'); c0=-1; T_l=2000; T_u=3000; N=50; end  
    if jj==36; data0=readmatrix('jobcorpY12.csv'); c0=-1; T_l=2000; T_u=3000; N=50; end      
    if jj==37; data0=readmatrix('jobcorpY1.csv'); c0=-1; T_l=1000; T_u=2000; N=50; end  
    if jj==38; data0=readmatrix('jobcorpY2.csv'); c0=-1; T_l=1000; T_u=2000; N=50; end      
    if jj==39; data0=readmatrix('jobcorpY3.csv'); c0=-1; T_l=1000; T_u=2000; N=50; end  
    if jj==40; data0=readmatrix('jobcorpY12.csv'); c0=-1; T_l=1000; T_u=2000; N=50; end      
    if jj==41; data0=readmatrix('jobcorpY1.csv'); c0=-1; T_l=40; T_u=1000; N=50; end  
    if jj==42; data0=readmatrix('jobcorpY2.csv'); c0=-1; T_l=40; T_u=1000; N=50; end      
    if jj==43; data0=readmatrix('jobcorpY3.csv'); c0=-1; T_l=40; T_u=1000; N=50; end  
    if jj==44; data0=readmatrix('jobcorpY12.csv'); c0=-1; T_l=40; T_u=1000; N=50; end  
    if jj==45; data0=readmatrix('jobcorpY1.csv'); c0=-1; T_l=40; T_u=3000; N=50; end  
    if jj==46; data0=readmatrix('jobcorpY2.csv'); c0=-1; T_l=40; T_u=3000; N=50; end      
    if jj==47; data0=readmatrix('jobcorpY3.csv'); c0=-1; T_l=40; T_u=3000; N=50; end  
    if jj==48; data0=readmatrix('jobcorpY12.csv'); c0=-1; T_l=40; T_u=3000; N=50; end      
    if jj==49; data0=readmatrix('jobcorpY1.csv'); c0=-1; T_l=2000; T_u=3000; N=40; end  
    if jj==50; data0=readmatrix('jobcorpY2.csv'); c0=-1; T_l=2000; T_u=3000; N=40; end      
    if jj==51; data0=readmatrix('jobcorpY3.csv'); c0=-1; T_l=2000; T_u=3000; N=40; end  
    if jj==52; data0=readmatrix('jobcorpY12.csv'); c0=-1; T_l=2000; T_u=3000; N=40; end      
    if jj==53; data0=readmatrix('jobcorpY1.csv'); c0=-1; T_l=1000; T_u=2000; N=40; end  
    if jj==54; data0=readmatrix('jobcorpY2.csv'); c0=-1; T_l=1000; T_u=2000; N=40; end      
    if jj==55; data0=readmatrix('jobcorpY3.csv'); c0=-1; T_l=1000; T_u=2000; N=40; end  
    if jj==56; data0=readmatrix('jobcorpY12.csv'); c0=-1; T_l=1000; T_u=2000; N=40; end      
    if jj==57; data0=readmatrix('jobcorpY1.csv'); c0=-1; T_l=40; T_u=1000; N=40; end  
    if jj==58; data0=readmatrix('jobcorpY2.csv'); c0=-1; T_l=40; T_u=1000; N=40; end      
    if jj==59; data0=readmatrix('jobcorpY3.csv'); c0=-1; T_l=40; T_u=1000; N=40; end  
    if jj==60; data0=readmatrix('jobcorpY12.csv'); c0=-1; T_l=40; T_u=1000; N=40; end  
    if jj==61; data0=readmatrix('jobcorpY1.csv'); c0=-1; T_l=40; T_u=3000; N=40; end  
    if jj==62; data0=readmatrix('jobcorpY2.csv'); c0=-1; T_l=40; T_u=3000; N=40; end      
    if jj==63; data0=readmatrix('jobcorpY3.csv'); c0=-1; T_l=40; T_u=3000; N=40; end  
    if jj==64; data0=readmatrix('jobcorpY12.csv'); c0=-1; T_l=40; T_u=3000; N=40; end        

if isnan(data0(1,1)); data0=data0(2:end,:); end       
T0=data0(:,3);    
idx_Ti=((T0>=T_l).*(T0<=T_u))>0;
data=data0(idx_Ti,:);
Y=c0*data(:,2);
X0=data(:,4:end);  
n=length(Y);
q1=round(n/N);
T1=data(:,3);  
Tmax=max(T1);
Tmin=min(T1);
T=(T1-Tmin)./(Tmax-Tmin);  % Transform T to unit interval
if TXs==0; X=X0; TX=[T,X]; p=size(X,2); end
if TXs==1
    X1=X0; X2=X0(:,max(X0)>1);
    X=[X1,X2.^2]; 
    TX=[T,X];
    p=size(X,2);
end

% nuisance parameter  
sn=sqrt(n);
an=0.15*log(n);
bn=0.85*(log(n)/(log(log(n))));
M=round(n^(2/3));

% Weight function Q
ln=0; for q=2:q1; ln=ln+q*(q-1)/2; end  
ind_ln=zeros(ln,4);
k=1; 
for q=2:q1
    for i=1:q
        for j=i+1:q   
            ind=[q,q*(q-1)/2,j,i];
            ind_ln(k,:)=ind;
            k=k+1;            
        end
    end
end
ind_q=(ind_ln(:,1).^(-2))./ind_ln(:,2);

% Cross-fitting
nK=floor(n/K);
Ki=1:1:K;
idx_Kn0=[kron(Ki',ones(nK,1));Ki'];
idx_Kn=idx_Kn0(1:n);
phat=zeros(n,1);
gamma=zeros(n,1);
B0_Ks=zeros(K,1);
B1_Ks=zeros(size(TX,2),K);

for k=1:K
    idx_k=(idx_Kn~=k);
    Y_k=Y(idx_k);
    T_k=T(idx_k);
    X_k=X(idx_k,:);
    n1=length(Y_k);
    idx_k0=(idx_Kn==k);
    T_k0=T(idx_k0);
    X_k0=X(idx_k0,:);
    
    % Conditional density esitmator p(T,X) 
    h_1=1.06*std(T_k)*n1^(-1/5); 
    lambda_global=1.1*norminv(1-(1/log(n1))/max(n1*h_1,p))*sqrt(n1);
    nK1=length(T_k0);
    phat0=zeros(nK1,nK1);
    parfor i=1:nK1
        warning('off')
        temp1=func_LogitCross(T_k<=(T_k0(i)+h_1),X_k,X_k0,n1,p,lambda_global,K_iter);
        temp2=func_LogitCross(T_k<=(T_k0(i)-h_1),X_k,X_k0,n1,p,lambda_global,K_iter);
        phat0(:,i)=max(trim,abs(temp1-temp2)/(2*h_1));
    end
    phat(idx_k0)=diag(phat0);     
    
    % Estimate gamma(t,x) by LASSO
    TX_k=TX(idx_k,:);
    TX_k0=TX(idx_k0,:);
    [B_lasso,FitInfo]=lasso(TX_k,Y_k,'CV',10);
    idx_cv=FitInfo.IndexMinMSE;
    B0_cv=FitInfo.Intercept(idx_cv);
    B1_cv=B_lasso(:,idx_cv);
    gamma(idx_k0)=B0_cv+TX_k0*B1_cv;
    B0_Ks(k)=B0_cv;
    B1_Ks(:,k)=B1_cv;
end

% Test statistics
idx=1;
zero_ln=zeros(ln,1);
Tn1_ln=zero_ln;
Tn2_ln=zero_ln;
vn_ln=zero_ln;
sig_ln=zero_ln;
psi_ln=zero_ln;
phi_ln=zeros(n,ln);

for q=2:q1
    % Instrumental functions
    r=1/q;
    qi=(0:r:1);
    qn=length(qi);
    ind_T1=(T*ones(1,q)<=(ones(n,1)*qi(2:qn)));
    ind_T2=(T*ones(1,q)>=(ones(n,1)*qi(1:q)));
    ind_T=(ind_T1.*ind_T2)>0;

    % numerical intergration of gamma(t,x)
    t_list=(0:1/(q*(M-1)):1)';
    tn=length(t_list);
    ind_t0=(M-1)*(0:q)+1;
    ind_t1=(t_list*ones(1,q)<=(ones(tn,1)*t_list(ind_t0(2:qn))'));
    ind_t2=(t_list*ones(1,q)>=(ones(tn,1)*t_list(ind_t0(1:q))'));
    ind_t=(ind_t1.*ind_t2)/M;
    gamma_int0=zeros(n,tn);
    for t=1:tn
        T_t=t_list(t)*ones(n,1);  
        TX_t=[T_t,X]; 
        for k=1:K
            idx_k0=(idx_Kn==k);
            TX_tk0=TX_t(idx_k0,:);
            gamma_int0(idx_k0,t)=B0_Ks(k)+TX_tk0*B1_Ks(:,k);
        end    
    end
    gamma_int=gamma_int0*ind_t;
  
    % Estimation of v(l) and sigma
    phi0=(((Y-gamma)./phat)*ones(1,q)).*ind_T;
    vi0=phi0+gamma_int;
    vn0=mean(vi0);    
    for i=1:q
        vn2=vn0(i);
        phi2=phi0(:,i)+gamma_int(:,i)-vn2;
        for j=i+1:q
            vn1=vn0(j);
            phi1=phi0(:,j)+gamma_int(:,j)-vn1;
            vn=vn2-vn1;
            phi=phi2-phi1;
            sig2=mean(phi.^2);
            if q==2 && i==1 && j==2
                sig0=sqrt(sig2);
                sig=sig0;
            else
                sig=max([sqrt(sig2),epsi*sig0]);
            end
            Tn1=sn*vn/sig;
            Tn2=max([Tn1,0])^2;
            psi=-bn*(Tn1<-an);  
            Tn1_ln(idx,1)=Tn1;
            Tn2_ln(idx,1)=Tn2;
            vn_ln(idx,1)=vn;
            sig_ln(idx,1)=sig;
            psi_ln(idx,1)=psi;
            phi_ln(:,idx)=phi;
            idx=idx+1;
        end
    end
end
Tn_KS=max(Tn1_ln);
Tn_CvM=Tn2_ln'*ind_q;

% Multiplier null distribution
cvB_KS=zeros(B,1);
cvB_CvM=zeros(B,1);
for b=1:B
    U=2*sqrt(3)*rand(n,1)-sqrt(3); 
    phi_u=(1/sn)*phi_ln'*U;     
    Tn1_u=phi_u./sig_ln+psi_ln;
    Tn2_u=max([Tn1_u,zero_ln],[],2).^2;    
    cvB_KS(b)=max(Tn1_u);
    cvB_CvM(b)=Tn2_u'*ind_q;
end
cvB_KS=sort(cvB_KS);
cvB_CvM=sort(cvB_CvM);
cvs_KS=quantile(cvB_KS,1-alphas+eta)+eta;
cvs_CvM=quantile(cvB_CvM,1-alphas+eta)+eta;
p_KS=mean(cvs_KS>=Tn_KS);
p_CvM=mean(cvs_CvM>=Tn_CvM);

%disp(['KS test: ',num2str(Tn_KS),'   ',num2str(p_KS)]);
disp(['CvM test: ',num2str(Tn_CvM),'   ',num2str(p_CvM)]);

clear gamma_int0 gamma_int vi0 phi_ln phi_u psi_ln sig_ln; 
clear Tn1_ln Tn1_u Tn2_ln Tn2_u vn_ln zero_ln phat0 phi0;
clear t_list ind_T ind_T1 ind_T2 ind_t ind_t1 ind_t2 ind_ln ind_q; 
clear data data0 TX_k TX_k0 TX_t TX_tk0 X0 X1 X2 X_k X_k0 Y_k;

time=toc;
disp(['time:    ',num2str(time)]);    

if jj==1; save job1108_Y01N50_T20T30P; end  
if jj==2; save job1108_Y02N50_T20T30P; end  
if jj==3; save job1108_Y03N50_T20T30P; end  
if jj==4; save job1108_Y12N50_T20T30P; end  
if jj==5; save job1108_Y01N50_T10T20P; end  
if jj==6; save job1108_Y02N50_T10T20P; end  
if jj==7; save job1108_Y03N50_T10T20P; end  
if jj==8; save job1108_Y12N50_T10T20P; end  
if jj==9; save job1108_Y01N50_T04T10P; end  
if jj==10; save job1108_Y02N50_T04T10P; end  
if jj==11; save job1108_Y03N50_T04T10P; end  
if jj==12; save job1108_Y12N50_T04T10P; end  
if jj==13; save job1108_Y01N50_T04T30P; end  
if jj==14; save job1108_Y02N50_T04T30P; end  
if jj==15; save job1108_Y03N50_T04T30P; end  
if jj==16; save job1108_Y12N50_T04T30P; end  
if jj==17; save job1108_Y01N40_T20T30P; end  
if jj==18; save job1108_Y02N40_T20T30P; end  
if jj==19; save job1108_Y03N40_T20T30P; end  
if jj==20; save job1108_Y12N40_T20T30P; end  
if jj==21; save job1108_Y01N40_T10T20P; end  
if jj==22; save job1108_Y02N40_T10T20P; end  
if jj==23; save job1108_Y03N40_T10T20P; end  
if jj==24; save job1108_Y12N40_T10T20P; end  
if jj==25; save job1108_Y01N40_T04T10P; end  
if jj==26; save job1108_Y02N40_T04T10P; end  
if jj==27; save job1108_Y03N40_T04T10P; end  
if jj==28; save job1108_Y12N40_T04T10P; end  
if jj==29; save job1108_Y01N40_T04T30P; end  
if jj==30; save job1108_Y02N40_T04T30P; end  
if jj==31; save job1108_Y03N40_T04T30P; end  
if jj==32; save job1108_Y12N40_T04T30P; end  

if jj==33; save job1108_Y01N50_T20T30N; end  
if jj==34; save job1108_Y02N50_T20T30N; end  
if jj==35; save job1108_Y03N50_T20T30N; end  
if jj==36; save job1108_Y12N50_T20T30N; end  
if jj==37; save job1108_Y01N50_T10T20N; end  
if jj==38; save job1108_Y02N50_T10T20N; end  
if jj==39; save job1108_Y03N50_T10T20N; end  
if jj==40; save job1108_Y12N50_T10T20N; end  
if jj==41; save job1108_Y01N50_T04T10N; end  
if jj==42; save job1108_Y02N50_T04T10N; end  
if jj==43; save job1108_Y03N50_T04T10N; end  
if jj==44; save job1108_Y12N50_T04T10N; end  
if jj==45; save job1108_Y01N50_T04T30N; end  
if jj==46; save job1108_Y02N50_T04T30N; end  
if jj==47; save job1108_Y03N50_T04T30N; end  
if jj==48; save job1108_Y12N50_T04T30N; end  
if jj==49; save job1108_Y01N40_T20T30N; end  
if jj==50; save job1108_Y02N40_T20T30N; end  
if jj==51; save job1108_Y03N40_T20T30N; end  
if jj==52; save job1108_Y12N40_T20T30N; end  
if jj==53; save job1108_Y01N40_T10T20N; end  
if jj==54; save job1108_Y02N40_T10T20N; end  
if jj==55; save job1108_Y03N40_T10T20N; end  
if jj==56; save job1108_Y12N40_T10T20N; end  
if jj==57; save job1108_Y01N40_T04T10N; end  
if jj==58; save job1108_Y02N40_T04T10N; end  
if jj==59; save job1108_Y03N40_T04T10N; end  
if jj==60; save job1108_Y12N40_T04T10N; end  
if jj==61; save job1108_Y01N40_T04T30N; end  
if jj==62; save job1108_Y02N40_T04T30N; end  
if jj==63; save job1108_Y03N40_T04T30N; end  
if jj==64; save job1108_Y12N40_T04T30N; end  

end

